Diversity-induced synchronized oscillations in close-to-threshold excitable elements 
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The question of how network topology influences emergent synchronized oscillations in excitable 
media is addressed. Coupled van der Pol-FitzHugh-Nagumo elements arranged either on regular 
rings or on clusters of the square lattice are investigated. Clustered and declustered rings are 
constructed to have the same number of next-nearest-neighbors (four) and a number of links twice 
that of nodes. The systems are chosen to be close-to-threshold, allowing global oscillations to be 
triggered by a weak diversity among the constituents that, by themselves, would be non-oscillating. 
The results clearly illustrate the crucial role played by network topology. In particular we found that 
network performance (oscillatory behavior and synchronization) is mainly determined by the network 
average path length and by the standard deviation of path lengths. The shorter the average path 
length and the smaller the standard deviation, the better the network performance. Local properties, 
as characterized by the clustering coefficient, are less important. In addition we comment on the 
mechanisms that sustain synchronized oscillations and on the transient times needed to reach these 
stationary regimes. 
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I. INTRODUCTION 



During the last decade, a considerable attention has 
been drawn to identify the mechanisms that mayinduce 
global oscillatory behavior in excitable media 0, 0, 0, 
IE IE IE IE IE IE EE IHI ■ Excitable media are composed 
of coupled excitable elements characterized by quiescent 
states with very large responses (amplitude and relax- 
ation time) to stimuli larger than some threshold. This 
kind of nonlinear behavior has attracted a great inter- 
est in fields such as chemistry 0, ^ , and biology . 
Oscillations in excitable media can be originated by: i) 
pacemaking cells such as in the case of the heart |l3| . 

ii) individual elements that oscillate when isolated, and, 

iii) more recently it has been suggested that diversity 
amongst the elements and/or internal noise can produce 
oscillations that are emergent in the sense that they show 
up even when none of the elements oscillate when iso- 
lated. 

Motivated by the experimental observation of emer- 
gent synchronized oscillations in the islets of Langerhans 
|LJ,|l5| (insulin producing structures in mammalian pan- 
creas) several works have been devoted to offer a sound 
explanation for this oscillatory behavior |E IE IE H^- 
Cells in the islets of Langerhans are of the type known 
as (3 cells and do not oscillate when isolated 0. As 
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no pacemaking cells are found in the islets of Langer- 
hans, a different mechanism has to be identified. In a 
recent work, Cartwright Q proposed that emergent syn- 
chronized oscillations may occur in a system of coupled 
excitable cells or elements in which half of them was silent 
and the other half continuously active, when taken indi- 
vidually. Each element was described by means of a van 
der Pol-Fitzhugh- Nagumo oscillator More recently 
we argued ^(| that, if the system was close to thresh- 
old, oscillatory behavior could be triggered by just a few 
diverse elements. In addition, for large enough coupling 
between the network elements, oscillations may become 
synchronized. This proposal is particularly interesting 
from a biological point of view as a weak diversity is in- 
herent to biological systems. 

In our previous work we arranged excitable elements on 
clusters of the square lattice. Such periodic two or three 
dimensional networks have been commonly used to de- 
scribe the topology of islets of Langerhans [E 0, 0] and 
some types of neurons [l9l | . The present work is addressed 
to investigate the role that network topology plays in the 
appearance of synchronized oscillations in an excitable 
medium described by the model proposed in [E . The 
effect of topology on synchronization of oscillators has 
been intensively investigated in recent years, particularly 
in networks having either the small-world or the scale- 
free property [IE HE EE EE Hi Most studies conclude 
that the larger the intrinsic size of the system (quantified 
by the average network length, see below) the worse the 
synchronization. In this work we investigate this issue on 
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regular networks with different average path lengths. We 
choose clusters of the square lattice with periodic bound- 
ary conditions and rings with additional links such that 
the degree k (number of nearest-neighbor connections) 
equals that of the square lattice clusters (fc=4). Our re- 
sults support the conclusion that the longer the average 
path length the worse the synchronization performance. 
In addition, they indicate that networks with a small 
standard deviation of the path length distribution show 
better synchronization. We also discuss the dynamical 
behavior of the coupled oscillators trying to identify the 
mechanisms that trigger and sustain synchronized oscil- 
lations. 

The paper is organized as follows. In Subscc. Ill Al 
we discuss the main features of the coupled van der Pol- 
FitzHugh-Nagumo equations outlining the main reasons 
why diversity may trigger oscillations. The various mea- 
sures we use to quantify the degree of oscillatory behav- 
ior and synchronization are described in Subsec. Ill Bl 
However, since we want to deal with (essentially) time- 
independent indicators, in the Appendix we discuss the 
stationarity of the present system which is dictated by 
the relaxation times to the limit cycles of the oscillators. 
In order to have a genuine indication of the stationarity, 
in the appendix we stick to standard indicators borrowed 
from circular statistics. The topologies of the networks 
investigated in this work are described in Subsec. Ill CI 
as well as the magnitudes used to characterize local and 
global structure, namely, the average path length and the 
clustering coefficient. The results are described and dis- 
cussed in Sec. IIIII We examine how oscillation and syn- 
chronization vary with the coupling between the elements 
in the network and the effects of the spatial distribution 
of diverse elements. The differences between networks 
are highlighted, aiming to identify the characteristics of 
network topology that determine network performance 
(as far as synchronization is concerned). We also con- 
sider the dependence of global oscillations on the size of 
the system. Finally, we address the question of how syn- 
chronized oscillations are sustained. The achievements of 
our research are summarized in Sec. IIVI 

II. MODEL AND METHODS 

A. Van der Pol-FitzHugh-Nagumo equations 

We describe the excitable medium by means of a sys- 
tem of coupled van der Pol-FitzHugh-Nagumo (PFN) 
equations |2^, |2(| 153, HI] , arranged on a given network. 
The set of equations is written as, 0, Il6l [2^ | : 

^ = 7 [zj - yf/3 + y t + kQ,] (la) 



-^ = -j- 1 (y l + v l + l3z l ). (lb) 



where i = 1, ,.,N, N being the total number of elements 
in the network. Like all relaxation oscillators, PFN oscil- 
lator has a variable yi with a fast release and a variable 
Zi with a slow accrual phase. The biological interpreta- 
tion of these equations is discussed in detail in |5j , here 
we briefly comment the main issues. Each oscillator is 
assumed to represent a cell in a living system. Variables 
yi and z% correspond to the potential across the nonlin- 
ear resistance (cell membrane) and the current through 
the supply, respectively. Parameters (3 and 7 describe the 
(membrane) resistance and the square root of the ratio in- 
ductance/capacitance respectively, and v is proportional 
to the charging potential of the cellular membrane. 

The strength of coupling between nearest-neighbors is 
denoted by k, while the function that describes coupling 
Qi is given by, 

fe 4 N 

Qt = J2(Vj -Vi) = -^2 Li i y i ' ( 2 ) 
i=i 3=1 

where ki is the number of nearest-neighbors (or connec- 
tivity) of element i. Lij is usually referred to as the 
Laplacian matrix, whose structure plays an essential role 
in determining the behavior of the network, as it contains 
all the information on network's topology; for instance it 
has been used to investigate the stability of the synchro- 
nized state once established |22j, |24j . 

In the absence of coupling (isolated element) the pa- 
rameter v determines the behavior of the oscillator. If v is 
outside the oscillating range (\v\ > O, O being the thresh- 
old), the element is in stable equilibrium. It is called an 
excitable element and can be either silent [y < — O) or 
continuously active (y > 0). Otherwise, the equilibrium 
is unstable (\v\ < 0) and the element performs oscilla- 
tions along the limit cycle. The threshold is given by: 

e = V7 2 - P ^3 (3) 

and for the chosen set of parameters 7 = 2 and (3 = 0.5, 
kept constant hereafter in this work, it takes the value 
= 0.60412. 

The coupled PFN equations were solved numerically 
on several regular networks with N nodes, 2N links and 
a connectivity fc, = 4 for all nodes. All elements are as- 
sumed to be active with v\ = v, but a small fraction x of 
them, called impurities or diverse elements, are taken to 
be silent with i>{ = —u, i = 1, . . . , [xN]. The numerical 
results discussed in this work correspond to v = 0.61. 
Although for this parameter value no element oscillates 
when isolated, coupling and impurities may trigger oscil- 
lations. A simple way to visualize this effect is to write 
(within a mean field approach) the effective parameter 
as: 



u eS = i/(l -2a:). (4) 
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This clearly shows that the presence of impurities may 
shift the effective parameter value below the threshold, 
forcing an originally quiescent element to oscillate. Al- 
though this expression can be used to estimate the frac- 
tion of impurities above which the oscillatory behavior 
emerges, it cannot account for the important differences 
in the performance of ordinary and impurity elements 
that can actually show up (see below). Fortunately, an 
effective v can be estimated separately for each element 
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At quiescence there is no coupling between elements hav- 
ing the same parameter value. However, in the case of 
one silent (uj — —v) and one continuously active cell 
(yi = +u) we can write (yj — t/j) « —(fj — h) = 2za 
The coupling between different cells is therefore of crucial 
importance for the effective change of cells parameters. 
Let us call m±(x) the average number of the nearest- 
neighbors of opposite nature for ordinary (+) and im- 
purity (— ) cells. The effective parameter can then be 
written as: 
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±v[l-2m±(x)0K] 
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For large enough N, the average number of nearest neigh- 
bors m± (x) is given by, 



m+(x) — n ( 4 _ x ) 

n=l ^ 
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B. Measures of oscillation and synchronization 

We quantify the emergence of oscillatory behavior by 
calculating various averages and standard deviations, 
over a sufficiently long time interval (see below), of ei- 
ther the fast variable yi or the phase, both for individ- 
ual elements and for the network as a whole. In the 
literature one can find many definitions of synchroniza- 
tion, mainly because of the different physical meanings 
(full-, phase-, delayed-synchronization, etc.) Measures 
similar to the ones adopted here are commonly used 
to characterize oscillatory behavior in excitable media 

HUHmEmiHHa. In order to make the brid s e 

with well-established circular statistics j^, in the Ap- 
pendix we evaluate the circular variance for some selected 
cases, drawing the same conclusions (at least qualita- 
tively) of this section as far as the level of synchronization 
is concerned. As a byproduct we also obtain information 
on the stationarity of our model. 

As discussed elsewhere 16] transients can be quite long 
for weakly cells. Instead, in the case of strongly coupled 
cells the coherent oscillatory behavior is almost instan- 
taneously attained (cp. Fig. 1 in Ref. [ly, see also 
Appendix). In order to minimize the contribution of 
transients [l6j|. the averages are computed starting from 
a large initial time t\ after which the stationary state 
is reached (for moderate and strong coupling regimes), 
and over a time interval At = tf — i; long enough to 
cover a sufficient number of periods. In a typical case, 
the coupled PFN equations were solved with an integra- 
tion step of 0.05, taking averages from t\ = 200 up to 
U = 400. The initial time t\ = 200 is large enough to 
surpass the transient period of moderately coupled net- 
works with low synchronization (see Appendix). For the 
values of cell parameters given above, the intrinsic period 
is around T = 10, so that oscillations were averaged over 
approximately 20 periods. 

We use the following standard deviation as a measure 
of intensity of cell oscillations (CO) at node j: 



This indicates that, on average, impurities can be much 
more affected by coupling than ordinary cells. To get a 
grasp of the implications of this result, let us give some 
numbers. The critical fraction of impurities x c required 
to force ordinary cells to oscillate can be derived by set- 
ting i/L ff = 9: 
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(9) 



Taking a typical value for the coupling strength k = 1 
and the values of the other two parameters given above, 
we obtain x c = 0.00241. The effective parameter of the 
impurity cell for this value of x is ~ +1.8. The effect 
is so strong that the parameter has not only changed its 
sign, but it is also well outside the oscillatory regime. 
This result underlies the differences in the behavior of 
ordinary and impurity cells discussed below. 
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t=ij 



where (yj) — 1/At^*f Vj(t) is the time average over in- 
terval At. A high value of crp° implies a large amplitude 
of cell oscillations around its time average, while a low 
value reveals an almost non-oscillatory behavior. 

We also define the average of cell oscillations over the 
whole network as: 
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Depending on the neighborhood, the effective parame- 
ters Vj S can largely differ from cell to cell (see preceding 
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subsec), inducing different oscillation intensities. The 
dispersion of cells oscillation intensity with respect to its 
network average can be quantified by means of the stan- 
dard deviation: 
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Synchronization is evaluated by comparing the current 
phases of each oscillator with the current phase of a ran- 
domly chosen reference cell (ordinary or impurity), and 
calculating the time average over the interval At [l6| : 



N t f 



aS = \ (N ~Tw ££ - cosMt)] 2 - (13) 

Shifting all the limit cycles to have a common center we 
have defined 



cos&(t) - m)/ri{t) = fe(i)/Vw(*) 2 + 5 i(*) a ( 14 ) 

where y } (t) = y 3 -(t) - (yj(t)) and Zj(t) = zj(t) - (zj(t)). 
The sum in <r s is not extended to all possible reference 
sites in order to limit the computational time. However, 
singling out a reference site allows us to track the dif- 
ference in the behavior of ordinary and diverse elements 
while, for instance, this information is hidden in a mea- 
sure like the circular variance that makes reference to an 
average phase (see Appendix) . The smaller a s the better 
the synchronization. The measure chosen to test synchro- 
nization in this work is more demanding than those used 
in previous works, in the sense that now we evaluate the 
differences in both dynamical variables 0, El • Note that 
the average values are subtracted from each dynamical 
variable in in Eq. (|13[) . This means that this measure 
does not account for possible differences in these aver- 
ages, that can arise from different effective parameters 
i/ off . 

In order to improve synchronization, first the major- 
ity of cells should be oscillating with the same frequency. 
Once that is achieved, dephasing between cells should 
vanish. Thus, an essential parameter for characterizing 
the emergence of synchronization is the dispersion of os- 
cillating periods. We quantify this dispersion by means 
of the standard deviation divided by the average period: 
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where T = l/iVEj=i 

is averaged over the whole 

network. 

Finally, we can measure the intensity of oscillations 
the network as a whole comparing the time dependent 



network oscillations with the temporal average over the 
interval At plj . 



NO 



\ 



(16) 



with the network activity y(t) = l/iVy^^ yj(t) and its 
temporal average (y) = l/At^t- y(t)- I n order to have 
large values of <r NO , the majority of cells must both os- 
cillate and be synchronized. Thus, the intensity of the 
network oscillations er NO is an appropriate quantity for 
the simultaneous measure of both activation and syn- 
chronization in excitable media 21:1 . 



C. Networks 

As the main aim of this work is to investigate the role 
of topology, we shall compare the results obtained on var- 
ious regular networks having a number of links twice that 
of nodes and the same degree (number of nearest neigh- 
bors). We choose the square lattice with periodic bound- 
ary conditions (this guarantees that all nodes i have de- 
gree ki — 4) and two types of regular rings (see Fig. 

The first one is the simplest clustered regular lattice, 
which consists of a ring with additional connections to 
the two next-nearest neighbors (this is usually referred 
as K = 2). The second type of a ring has zero clustering 
(see below) and is constructed as follows. Standard (clus- 
tered) rings with K — n are built by linking each site to 
all neighbors from the first up to the n-th. Declustered 
rings (see |3^|), denoted as K = n, have additional links 
only to the n-th neighbors. In this work we investigate 
rings with K = 2 and K = 3, both having degree fc^=4 
for all nodes. 

The topology of connections between sites can strongly 
influence the cooperation and the exchange of informa- 
tion between them. The structural properties of a graph 
are usually quantified by the average path length L and 
the clustering coefficient C 35, 36|,|37|]. The average path 
length is calculated as the network average of the shortest 
graph distances between two nodes for all possible pairs, 



L = 



N(N-l)^ dlJ - 



(17) 



The clustering coefficient C measures to which extent 
are the neighbors of each site connected to each other. It 
is defined as follows, 



C 



AT X/X/ 



(18) 



where j and k are nearest-neighbors of i and rijk = 1 only 
if there is a link betwee j and k. This definition does only 
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FIG. 1: Illustrates the link structure in clustered rings (upper) 
with connections to the first and second nearest-neighbors 
(K — 2), and declustered rings (lower) with connections to 
the first and third nearest-neighbors (K = 3). Note that while 
the clustered ring has triangles of elements the declustered 
ring has quadrilaterals. 



give C 7^ if there are triangles in the network. This 
is the case of clustered rings such as that chosen here 
(K = 2) but it is not that of the square lattice nor that 
of the declustered ring, both having C = 0. 

The average length for the three networks investigated 
here can be calculated analytically. The result for the 
clustered ring is, 



L(K 



N N/2 + 1 
T N- 1 



(19) 



valid only for N — An, where n is an integer. For very 
large systems this behaves as L ~ N/8. In the case of 
the declustered ring we obtain, 



L(K = 3) 



N N/2 + A 
~6 N- 1 



(20) 



valid only for networks with a number of nodes N — 6n, 
n being an integer. In this case the limit for N — ► oo is 
L ~ N/12. 

It is interesting to write a general formula for an arbi- 
trary K in the case of a declustered ring, 



L(K) = 



N N/2+ (K 2 -i)/2 



2K 



N -1 



(21) 



where £=1,2 for K odd or even, respectively. This for- 
mula is valid for N = 2Kn, where n is an integer. 

Finally, the result for clusters of the square lattice with 
periodic boundary conditions is, 



L (square) 



AT3/2 



2(N - 1) 



(22) 



its asymptotic behavior being L ~ 0.5y/N. 

The main difference between rings (clustered or declus- 
tered) and the square lattice is that the latter is far 
smaller than the former as its average length increases 
sublincarly with the network size. It is worth trying 
to understand why the square lattice behaves in this 
way in terms of the structure of the Laplacian matrix. 
By appropriately renumbering sites, any cluster of the 
square lattice with periodic boundary conditions can be 
mapped onto a ring in which additional links are intro- 
duced. These links are: links between y/N nodes and 
their (y/N — l)-th and (y/N + l)-th neighbors, and links 
between the remaining N — y/N and their (y/N + l)-th 
neighbors. For large N this is equivalent to a declustered 
ring with K = y/N . In fact Eq. I|2U|) for declustered rings 

gives, for K = y/N, the same assymptotic behavior than 
that obtained for clusters of the square lattice (see Eq. 
(12211 1. Contrary to regular rings with K — 2 or K = 3 
where additional links introduce non-zero elements in the 
second and third subdiagonals of the Laplacian matrix, 
respectively, links in the square lattice are long-ranged 
and introduce non-zero elements in the (y/N — l)-th and 
(y/N + l)-th subdiagonals. 



III. RESULTS 

Results discussed hereafter were obtained for the fol- 
lowing values of the model parameters: 7 = 2 and j3 = 0.5 
were kept constant, while we took v = 0.61 for normal 
elements and v = —0.61 for diverse elements. The con- 
centration of impurities is 2%. The number of elements 
was mainly fixed at N — 100 although we also consider 
clusters of the square lattice of size 18 x 18. In all cases 
initial conditions were yi(0) = 1 and 2^(0) = 1 for all net- 
work nodes i. We have previously checked that varying 
the initial data within the limit cycle has a very small 
effect on the indicators of Subsec. IIIBI evaluated after 
the transients. 



A. Coupling Strength 

The analysis of the influence of coupling between cells 
on intensity of oscillations and synchronization reveals 
three different regimes (see Fig. [2J. The results shown in 
the Figure correspond to a cluster of the square lattice 
having N = 100 nodes with a fraction of impurities x = 
0.02. 

In the case of weak coupling (k < 0.1), the chosen 
transient t\ = 200 is not long enough to allow the system 
to reach the stationary state. The behavior of the sys- 
tem is strongly influenced by the initial conditions. The 
improvement of synchronization as coupling strength de- 
creases to zero is a consequence of the initial conditions 
(common to all elements). Increasing t\ would allow the 
weakly coupled oscillators to run away from each other, 
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FIG. 2: Measures used to characterize the emergence of os- 
cillations and synchronization versus strength of coupling be- 
tween excitable elements in the PFN system described by Eqs. 
Q. Excitable elements are arranged on a square network of 
linear size I = 10, with 98 continuously active sites (y — 0.61) 
and 2 silent impurities (y = —0.61). The results are aver- 
aged over 10 different realizations of impurity distributions. 
Cells oscillation (cr co , see Eq. Ull^ - ) - thick dot-dashed line, 
network oscillations (<7 NO of Eq. 1161 ') - dotted, and synchro- 
nization (<r s , see Eq. 1)13^ with respect to an ordinary (full), 
or an impurity element (dashed). Dispersion of element os- 
cillation (cr DO , of Eq. C2J) - thin dot-dashed, and of the 
oscillation periods (cr DP , see Eq. 1151 ') - stars (see text) are 
also shown. 



FIG. 3: Synchronization (circles) a s of Eq. 113H . dispersion of 
oscillation (squares) <j do of Eq. Ill 2ft , and dispersion of peri- 
ods (stars) cr DP of Eq. (JUJ, versus cell oscillation <r co of Eq. 
lllll . Elements were arranged on clusters of the square lattice 
with N — 100 with 98 continuously active sites [y = 0.61) and 
2 silent impurities (y — —0.61). The results are averaged over 
10 different realizations of impurity distributions. Coupling 
strength is varied over the range n £ (0, 1). 

only in the case of very strong coupling, cp. Fig. [5J Rea- 
sons for the different behavior of ordinary and impurity 
elements were given in Subsec. II. A. 



see Appendix, deteriorating synchronization. Only after 
rather long t\ the network will reach the stationary oscil- 
lating state. As k goes to zero the transient period will 
become infinite. Results shown in Fig. for such low 
values of coupling cannot be considered too reliable. 

At moderate coupling (0.1 < k < 0.3), the pertur- 
bation originating from impurity cells is still not spread 
much over the system. Ordinary cells are only weakly 
impelled to oscillate, but with different amplitudes. The 
dispersion of oscillations is large, reaching its maximum 
for k « 0.15, cp. Fig. [21 In addition, a moderate cou- 
pling is not strong enough to eliminate large dephasing 
between elements oscillating with the same frequency. 
Consequently, the oscillation of the network as a whole 
is much weaker than that of individual elements, as syn- 
chronization is still not attained. 

In the strong coupling regime (n > 0.3), oscillations 
of elements become more and more intensive and syn- 
chronization is gradually improved, leading to a strong 
oscillation of the network as a whole. Elements inten- 
sively oscillate and there are no significant differences in 
frequencies (the dispersion of periods is practically zero). 
Increasing the coupling strength further, the phase dif- 
ferences are eliminated and synchronization is improved. 
Finally, ordinary elements are much better synchronized 
with each other, while the impurity elements are forced 
to oscillate in the same way as the rest of the network 




coupling strength 

FIG. 4: Same as in Fig. |2]but with different network topolo- 
gies: Excitable elements are arranged either on clustered (up- 
per) or declustered (lower) rings both with iV=100 nodes, 
with 98 continuously active sites (y — 0.61) and 2 silent im- 
purities (y = —0.61). 

The aforementioned mechanism of improving synchro- 
nization is even more evident if we plot intensity of os- 
cillations, network synchronization, and dispersion of in- 
tensities and period of oscillations in a single graph (Fig. 
EJ. Synchronization is the worst when the deviation of 
periods reaches its maximum for er co around 0.3. After 
that, it is practically unchanged until all elements oscil- 
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late with the largest amplitude. Approaching the max- 
imal activation, the dispersion of oscillations decreases, 
and synchronization is rapidly improved by eliminating 
differences in phases. 

In the case of rings, the values of coupling strengths 
needed both to shorten the transient period below a given 
t{ and to synchronize the elements are larger, cp. Fig. 
0] The long average path length of the clustered ring 
(L = 12.88 to be compared with L = 5.05 for square 
lattice of the same size), hinders the effective spread of 
perturbations coming from impurity elements. Both the 
dispersion of oscillations and the dispersion of periods 
remain large. The network reaches the steady regime 
for k > 0.5. Similar results are obtained for declustered 
rings, with a modest improvement of er s . Oscillation and 
synchronization are again quite low, due to the large av- 
erage path length (L = 9.09), cp. Fig. 0] This indi- 
cates that collective oscillatory behavior is determined 
mainly by the average length of the network, and not 
by its local properties quantified by the clustering coeffi- 
cient. Note that the square network and the declustered 
ring, both having zero clustering and appeciably different 
average path lengths, have very different performances. 
The shorter the length the better the performance. 
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Fig. shows the probability density function of path 
lengths in these networks. The density shows a well de- 
fined peak at the average path length in the case of the 
square lattice cluster, while it is almost constant in the 
declustered ring. As a consequence the standard devia- 
tion of path lengths is appreciably smaller in the former 
(3.7 to be compared with 4.9 for the declustered ring). 
We believe that the different standard deviations of path 
lengths is the origin of the noticeably different perfor- 
mance of the two networks. This result suggest that more 
homogeneous systems, as far as path length is concerned, 
are more efficient in triggering collective oscillations and 
synchronization. 
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FIG. 6: Probability density function of path lengths in a 
declustered ring with iV=100 nodes and K = 3 (continuous 
line) and in a cluster of the square lattice with N = 18 2 nodes 
(broken line). 



FIG. 5: Same as in Fig. [5] but with a different amount of 
diversity: Excitable elements are arranged on a square net- 
work of linear size I = 18, with 318 continuously active sites 
(y = 0.61) and 6 silent impurities (y — —0.61). 

Now, the question is whether the average path length 
is or not the only topological feature that matters. We 
address this question analyzing a cluster of the square 
lattice of a bigger size (18x18, with N — 324), hav- 
ing almost the same characteristic length (L = 9.03) 
as a declustered ring with K = 3 and 100 sites (see 
above). The coupling strengths are varied over the same 
range (k <e (0, 1)), while the fraction of impurities is not 
changed (ir=0.02 as above). As shown in Fig. [5] per- 
formance is now much better than that found for the 
declustered ring (see Fig. 0J| indicating that the average 
length is not enough to characterize network behavior. 



B. Distance between diverse elements 

All the results presented up to now are averaged over 
several realizations of impurity distribution. It turns out 
that the oscillating behavior strongly depends on the dis- 
tance between impurity elements. It is likely that larger 
distance between them would improve the performance, 
as the overlapping of their respective influence domains 
is smaller. In other words, the average path length from 
ordinary sites to the closest impurities becomes shorter. 
This is illustrated in Fig. [7| where synchronization is 
plotted as a function of the distance between impurities. 
Rings and square clusters behave in a similar way (the 
pace at which synchronization improves with distance is 
similar in all networks) although synchronization is much 
poorer in the former. 
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FIG. 7: Synchronization (ordinary reference site) versus dis- 
tance between the two impurity elements in several networks 
of size TV = 100: square network (square), clustered (trian- 
gle) and declustered (diamond) rings. Coupling strength and 
fraction of impurities are k = 0.5 and x = 0.02, respectively. 

C. The Size of the Network 

We have investigated how synchronization varies with 
the network size in the case of the square lattice. Calcu- 
lations were carried out keeping the fraction of impurities 
constant x = 0.02. As argued in 8] the effective coupling 
constant varies with the average network length as 1/L 2 . 
This can be easily checked by noting that in the con- 
tinuum limit the coupling term of Eq. (J2J behaves as, 



FIG. 8: Element oscillation (dot-dashed line), network oscil- 
lation (dotted line) and synchronization with respect of an 
ordinary reference site (full line) versus size of clusters of the 
square lattice. The coupling strength was taken either con- 
stant (k = 0.5, thin lines) or proportional to the square of the 
average path length (k = 0.5./V/100, thick lines). 

constant, again due to the fast (linear) increase of the 
average length with the network size. 



D. Sustaining synchronized oscillations 



Q(X 5#' °<V<L. (23) 

Then, rescaling the spatial variable 77 to bring to evidence 
size effects we find, 

1 d 2 y 

Q(X ZW' °<V<1, (24) 

where L should be interpreted as the average linear size 
of the system. 

In order to check this behavior we have carried out 
calculations with a coupling strength either constant 
(k = 0.5) or proportional to L 2 (k = 0.5JV/100). The 
numerical results are depicted in Fig. [5] We first note 
that, for k = 0.5, synchronization is gradually deterio- 
rated as the size of the system increases. We have checked 
that the chosen transient t\ = 200 is long enough for the 
largest network shown in the Figure N = 400. The pace 
at which it worsens is rather slow due to the sublinear 
behavior of the average path length (see above) . Rescal- 
ing the coupling strength as just mentioned, gives a syn- 
chronization that remains almost unchanged as the size 
increases, proving the validity of the scaling argument. 
A similar analysis for rings shows a much faster wors- 
ening of the system performance when coupling is kept 
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FIG. 9: Time dependence of oscillation of a randomly cho- 
sen ordinary element (yi(t) - full line), its coupling function 
(10 x Qi(t) - dot-dashed) and oscillation averaged over the 
whole network ((y(t)) - dots almost superimposed to yi(t)). 
Elements were arranged on a square cluster of linear size 
I = 10. The results correspond to a single realization of im- 
purity distribution. The chosen coupling strength (k — 0.5) 
is large enough to trigger synchronized oscillations. 

The results and discussion of the preceding subsections 
clearly show that coupling between elements and diver- 
sity are the main cause of synchronized oscillations in the 
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systems under study. Coupling is finite only if yi are not 
all identical. Then, how is synchronization sustained? 
Fig. El shows yi(t) and the coupling function Qi(t) for a 
single element i as a function of time. The average over 
the whole network (y(t)) is also shown. The results cor- 
respond to the strong coupling regime (k — 0.5). First 
we note that a characteristic oscillating pattern and co- 
herent behavior is achieved almost instantaneously. We 
have checked that this is the case even if the initial con- 
ditions are varied. On the other hand no significant dif- 
ferences between oscillation of a particular element and 
of the network as a whole are noted. The coupling func- 
tion Qi (t) shows an oscillatory pattern similar to that of 
yi, although it is an order of magnitude smaller. This 
shows that a small difference between variables in each 
element is enough to sustain synchronized oscillations. 
In addition the oscillatory pattern is stable: any pertur- 
bation of phase coordinates (within the range of limit 
circle) at time £ , could be viewed as setting new initial 
conditions starting at that time. The effect of the per- 
turbation would then be eliminated after a few periods, 
restoring the characteristic pattern. 



IV. CONCLUSIONS 



creases, provided that the coupling strength is increased 
proportionally to the square of the average path length. 
Regarding synchronization, the optimal topologies are 
the ones for which the average path length increases more 
slowly with the size of the system. 
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APPENDIX A: CIRCULAR VARIANCE AND 
STATIONARITY 

The degree of synchronization of the coupled PFN os- 
cillators as well as its dependence on time can be anal yzed 
using well-established measures of circular statistics |34| . 
Given the set of phases 4>i(t) of Eq. l|14fl as a function of 
time and the corresponding unit vectors OPi, the average 
phase <5(i) is defined to be the direction of the resultant 
OP. The cartesian coordinates of P are: 



The conclusions that emerge from our study of the 
induction of global synchronized oscillations by hetero- 
geneity in excitable regular systems are: 

1. Elements and network activity and synchronization, 
as collective phenomena, depend mainly on global prop- 
erties of the network. Shortening the average path length 
by merely changing the topology, while keeping constant 
the number of sites and links, improves synchronization. 
In addition we found that the standard deviation of path 
lengths is also a very crucial variable: the smaller it is 
the better the network performance becomes. This in- 
dicates that homogeneity, as related to path length, is a 
benefitial feature of networks. 

2. Local properties quantified by the clustering coeffi- 
cient do play a role much less significant than the average 
path length and its standard deviation. 

3. Synchronization in clusters of the square lattice is 
improved by increasing the strength of coupling between 
the excitable elements. In the strong coupling regime 
most of the elements are highly activated, oscillating with 
the same frequencies and negligible dephasing. 

4. Circular regular networks with short range links 
are hardly synchronized even when elements are strongly 
coupled, due to their long average path length and large 
standard deviation. 

5. Oscillation and synchronization depend on the spa- 
tial distribution of impurity elements, improving as the 
distance between those elements increases. 

6. Finally, we have revisited the behavior of van der 
Pol-FitzHugh-Nagumo excitable media when their size 
is varied. We found that the dynamical behavior of the 
network can be kept unchanged when the system size in- 



C(t) = -Y, C0S Mt), S(t) = -5>in<Mi), (Al) 

i=l i=l 

leading to the length of the resultant: 

R{t) = NR(t) = Ny/C(t) 2 + S{t) 2 . (A2) 

The average phase $(t), i.e. the mean direction, is 
obtained from the following equations: 

C{t) = R{t) cos $(t), S(t) = R(t) sin $(t). (A3) 

Recalculating the current phases relative to the mean 
(£(t) = 'Pit) — &(t)), one can see that the sum of devia- 
tions about the mean equals zero: 

N N 

5>m&(*) = $>in(&(i) - *(<)) = 0. (A4) 

i=l i=l 

As a possible measure of the synchronization one could 
use the circular variance [34j], given by: 

1 - 

tf(*)=jy£[l-«»(&(t)-<lKi))]. (A5) 

1=1 

It can be shown that the variance D is minimized, if the 
phases are recalculated relative to the mean |3^| . 

Small values of D indicate that the phases of coupled 
oscillators are grouped around the mean, hence oscilla- 
tors being well synchronized. On the other hand, large 
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FIG. 10: Time dependent circular variance D(t) for N = 100 
oscillators with two impurities coupled into a square lattice 
with strong and moderate interaction. Upper: re = 0.5; 
Lower: re = 0.125. 
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FIG. 11: Time dependent circular variance D(t) for N = 100 
oscillators with two impurities coupled into a square lat- 
tice with small interactions. Upper: re = 0.025; Lower: 
re = 0.0125. 



values of D correspond to rather uniform distribution of 
phases and low synchronization. 

Let us analyze first the square lattice. In the case of 
strong coupling (k = 0.5) oscillators are well synchro- 
nized, with D(t) having very small values, cp. Fig. ITUI 
(upper). Better synchronization is achieved in the slower 
part of the limit cycle, while oscillators tend to run away 
from each other in the faster part of the path, increasing 
D(t) towards 0.045. The coupling is strong enough to 
lead to stationary behavior already after 2 or 3 periods. 

Decreasing the coupling strength, deteriorates syn- 
chronization. In the case of k = 0.125, the variance D(t) 
takes much larger values (up to 0.5 in stationary regime), 
while the transient extends to t = 125. Decreasing fur- 




time 

FIG. 12: Time dependent circular variance D(t) for N = 100 
oscillators coupled into a regular ring with two impurities and 
re = 0.5. Upper: clustered ring, K = 2; Lower: declustered 
ring, K = 3. 

ther the coupling strength (k = 0.025) worsens the syn- 
chronization and increases the transient period (t=500), 
cp. Fig. 1111 flower') . In the case of ultra- weak coupling 
(re = 0.0125), the oscillators are at the beginning quite 
well synchronized (for t < 300), as their behavior is still 
dominated by the common initial conditions. Later on, 
the phase differences are gradually increased, leading to 
almost uniform distribution for D(t) > 0.9. 
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FIG. 13: Time dependent circular variance D(t) for N = 324 
coupled oscillators arranged on a cluster of the square lattice 
with six impurities and re = 0.5. 

In the case of regular rings, better synchronization 
is achieved for declustered rings, having shorter aver- 
age path length, cp. Fig. ^] The average variance is 
smaller, 0.4 for declustered comparing to 0.6 for clustered 
systems, while the transient period is shorter, t « 100 
for declustered relative to t « 200 for clustered rings. 
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As reasonable, from the numbers above we see that the 
transient time necessary to attain a stationary evolution 
increases when the coupling constant is decreased. 

Finally, as in the main text, we address the question of 
whether the average path length is the only topological 
parameter that matters. To this end, we compare the 
circular variance for the declustered regular ring of 100 
elements with that for a cluster of the square lattice with 
324 elements. Both networks have squares as basic cy- 



cles and almost the same average path lengths, while the 
standard deviation of the distribution of path lengths is 
smaller for the square lattice (see Subsec. MI Ajl . Com- 
paring Figs. 1131 and IT21 flower) we note that the square 
network has a smaller average variance (D w 0.2) and a 
shorter transient (t « 60). Thus, we conclude again that 
smaller deviation of path lengths improves synchroniza- 
tion and shortens the time needed to reach the steady 
state. 
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